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The spin decomposition of the proton is a long-standing topic of much interest in hadronic physics. 
Lattice QCD has had much success in calculating the connected contributions to the quark spin. 
However, complete calculations, which necessarily involve gluonic and strange-quark contributions, 
still present some challenges. These “disconnected” contributions typically involve small signals 
hidden against large statistical backgrounds and rely on computationally intensive stochastic tech¬ 
niques. In this work we demonstrate how a Feynman-Hellmann approach may be used to calculate 
such quantities, by measuring shifts in the proton energy arising from artificial modifications to 
the QCD action. We find a statistically significant non-zero result for the disconnected quark spin 
contribution to the proton of about —5% at a pion mass of 470 MeV. 


I. INTRODUCTION 

The simple quark model picture suggests that the to¬ 
tal nucleon spin is comprised entirely in terms of its con¬ 
stituent quark spins. In contrast, experimental measure¬ 
ments reveal that the quark spin only generates about 
one third of the total nucleon spin [1]. This observation 
is a clear representation of the nontrivial dynamics as¬ 
sociated with nonperturbative QCD. Resolving the full 
composition of the nucleon spin in terms of the QCD 
degrees of freedom remains an active experimental and 
theoretical pursuit. For an overview of the status and 
progress, we refer the reader to Refs. [IHZl- 

As a systematically improvable method for studying 
nonperturbative properties of QCD, lattice simulations 
offer the potential to provide quantitative predictions for 
the decomposition of the nucleon spin. For recent numer¬ 
ical investigations of the nucleon spin, and related matrix 
elements, see Refs. [HHISl- 

In the conventional approach, spin matrix elements 
are extracted from 3-point correlation functions. Oper¬ 
ator insertions that are directly connected to the quark 
field operators of the nucleon interpolators can be reliably 
computed using established techniques. The operator in¬ 
sertions that involve self-contracted fermion lines, which 
are essential to isolate the strangeness spin content, for 
instance, require the stochastic estimation of the trace of 
an all-to-all propagator. Owing to the increased compu¬ 
tational demand of this stochastic estimator and a rel¬ 
atively weak numerical signal, such disconnected contri¬ 
butions have often been neglected in lattice simulations. 
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Nevertheless, there has been substantial progress made 
in recent years UMI]. For a related calculation involv¬ 
ing the vector current matrix elements we also refer to 

Ref. dg. 

In recent work, we have proposed an alternative to the 
conventional 3-point function technique for the study of 
hadron matrix elements in lattice QCD. By adapting the 
Feynman-Hellmann (FH) theorem to the lattice frame¬ 
work, we are able to isolate matrix elements in terms of 
an energy shift in the presence of an appropriate weak 
external field [23l El]. This is similar to the technique 
proposed by Detmold in Ref. E^. In Ref. E3| we used 
the Feynman-Hellmann relation to extract the gluonic 
contribution to the nucleon mass. The application of 
Feynman-Hellmann was further developed in Ref. |24j for 
the study of the connected spin contributions in various 
hadrons. We have also recently shown how it is possi¬ 
ble to compute flavour-singlet renormalisation constants 
nonperturbatively by an appropriate application of the 
FH theorem [55]. 

In the present work, we apply the FH technique to re¬ 
solve disconnected spin matrix elements. Whereas the 
connected spin contributions could be computed on con¬ 
ventional gauge helds, the disconnected contributions re¬ 
quires the generation of new special-purpose gauge con¬ 
figurations. The influence of the weak external spin 
field is therefore accumulated through the importance 
sampling of the hybrid Monte Carlo simulation. While 
such new configurations come at significant computa¬ 
tional cost, the computing time is comparable to that 
required for reliable with sampling using the conventional 
stochastic techniques. 

The manuscript proceeds as follows: Section H reviews 
the implementation of the FH theorem for the extraction 
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of spin matrix elements and summarises the lattice sim¬ 
ulation parameters; Section III describes the strategy for 
the isolation of the relevant matrix elements from the 
two-point correlation functions; with results reported in 
Section IV; followed by concluding remarks in Section V. 


II. FEYNMAN-HELLMANN METHODS AND 
SIMULATION DETAILS 


Here we discuss the Feynman-Hellmann approach in 
the context of calculations of disconnected contributions 
to matrix elements, in particular the quark axial charges. 
For details of previous calculations of the connected con¬ 
tributions, and the Feynman-Hellmann technique in gen¬ 
eral, see [24] , 

The quark axial charges are defined by forward matrix 
elements of the axial operator, 

(p,s| 9 ( 0 ) 7 ^ 759 ( 0 ) Ip, s) = 2is^Aq. ( 1 ) 

We access disconnected contributions to these quantities 
by implementing a modification to the fermion part of the 
QCD Lagrangian during gauge-field generation. Extra 
terms are included involving the axial operator weighted 
by some freely-varying real parameter A, applied equally 
to all three quark flavours, 


C-)- C + X ^ 973759 • (2) 

q—u,d,s 

This operator satisfies 75 -hermiticity, and so the determi¬ 
nant of the fermion matrix is still real. Hence we avoid 
introducing any sign problems. We choose projection 
matrices to isolate spin-up and down components of the 
nucleon correlation function, 

r± = ^(I + 74)^(II±i7573), (3) 

and by application of the Feynman-Hellmann relation, 
find that the correlator picks up a complex phase which 
mimics an imaginary energy component, 

E E{X) + i(j){X). (4) 


At first order in the parameter A, there is no shift in 
the real part of the energy, and the shift in the phase 
is exactly equal to the disconnected contribution to the 
total quark axial charge. 


dE 


= 0 

A=0 


d(j) 

dX 


= iAEdisc. , 

A=0 


( 5 ) 


where the total contribution is the sum of the individual 
flavour contributions. 


ASfiisc. — Audisc. T Addisc. T As. (6) 


{ki,Ks) X 

-^conf. 

^sources 

<(> 

-0.0125 

500 

1 

0.0014(10) 

(0.120900, 0.120900) -0.00625 

500 

1 

0.00002(83) 

0.03 

500 

1 

-0.00237(77) 

(0.121095, 0.120512) 

600 

1 

-0.0008(13) 

0.05 

800 

5 

0.00027(61) 


TABLE I. Table of ensembles generated for this work. Two 
pion masses with three and two values of A respectively have 
been used. The number of configurations and sources used, 
as well as the phase shift measured (discussed in Sec. Ill and 
Sec. IVI are also listed. 


Also note that the strange contribution is purely discon¬ 
nected. The different signs in Eq. ([^ result from the 
different choices of r±, and we note that changing the 
spin projection is equivalent to flipping the sign of A. 

Our strategy for this calculation, motivated by Eq. ([^ , 
is to generate new gauge ensembles for multiple values 
of A, measure the phase shift in Eq. Q and determine 
AEdisc. from the linear behaviour. 

In our previous work, we were able to access the con¬ 
nected part by implementing the change in Eq. to the 
Dirac matrix before inversion to compute the quark prop¬ 
agator entering hadron correlation functions (see |24 ) ). 
Here the modification is made to the fermion matrix in 
the HMC algorithm, and so information about the purely 
disconnected contributions is accessed. 


A. Simulation details 

We use gauge field configurations with 2-1-1 flavours of 
non-perturbatively 0(a)-improved Wilson fermions and 
a lattice volume of xT = 32^ x 64. The lattice spacing 
a = 0.074(2) fm is set using a number of singlet quantities 
[27H30] . The clover action used comprises the tree-level 
Symanzik improved gluon action together with a stout 
smeared fermion action, modified for the implementation 
of the Feynman-Hellmann method |24j . 

For the results discussed here, we use ensembles 
with two sets of hopping parameters, {ki,Ks) = 
(0.120900,120900) and (0.121095,0.120512), correspond¬ 
ing to pion masses of approximately 470 and 310 MeV. 
These have been generated with the modified quark ac¬ 
tion described in Eq. ([^. The details of these ensembles, 
including the values of A realised, are given in Table |Tj 


III. ANALYSIS TECHNIQUES 

A standard zero-momentum projected nucleon corre¬ 
lation function is given by 


Note that we access the total contribution because the 
operator in Eq. ([^ includes terms for all 3 quark flavours. 


G± (t) = J d^x r± (n I N{x)N{0) I n ) -1^^^ Ae-®*, 
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where N and N are interpolating operators coupling to 
the nucleon ground state, and the projection matrices 
F-i- (defined in Eq. (|3|) project spin-up and down com¬ 
ponents respectively. For our simulations, we use identi¬ 
cal source and sink smearing and operators. Hence, the 
amplitude A is purely real. 

With the modification to the Lagrangian in Eq. ([^, 
an imaginary component is introduced to the exponential 
factor, in addition to a complex shift in the amplitude. 
This shift in the amplitude is not the focus of this work, 
but is related to the A dependence of the wavefunction 
overlap factors. To first order in A, the amplitude and 
energy take the form 

A + X{AA + iAB), (7) 

E^E + iXAX:, (8) 




t/a 


and the correlation function at large times is given by 


G±{X,t) [A±X{AA + iAB)]e-^^^^^^^^G (9) 

(Recall that changing the spin projection corresponds to 
flipping the sign of A, as discussed in Sec. 0 - The quan¬ 
tity of interest is the shift in the phase, AE. To extract 
this value, we form the following ratio of real and imagi¬ 
nary parts of spin-up and down projections, 


Im[G-(A,t)-G+(A,t)] 
Re[G_(A,t) + G+(A,t)] 
large* sin (A AEt) — A ^ COS (A AE<) 
cos (AAEt) -I- A^ sin (AAEt) 


( 10 ) 


Note that the form of this ratio does not change if we in¬ 
clude second order terms in Eq. 0 and Eq. ([^. For the 
operator included in Eq. ([^, we find that the second or¬ 
der shift in the energy is purely real, and it can be shown 
that only the factor will change. Hence, corrections 
to these calculations do not appear until O (A^). 

The ratio in Eq. (101 is what we fit in our analysis. To 
determine ground state saturation of this quantity, we 
observe that, provided t ttie behaviour of the 

ratio is approximately linear in t. 


1 

R{X, t) = AAEt — A—^ -1- O (A^) , a |AA^ ' 

( 11 ) 

Previous determinations of AE [IMI] suggest that we 
should expect |AE| « 0.1, and hence for the largest value 
of A realised on our ensembles, aX = 0.05, this linear 
approximation will hold for times | <C 200. With this, 
we are able to introduce an ‘effective phase shift’. 


= ^[7^(A,^ + a)-7^(A,^)] (12) 

which in the regime discussed has the behaviour 
(jieff. = AAE , a < t < 


FIG. 1. Plots of the ratio in Eq. ( |10[ ) and the effective phase 
shift defined in Eq. (131 for A = 0.03, m,r ~ 470 MeV. The 
fitting window (shown in darker blue) was between time slices 
2 and 12. The errors shown are from a bootstrap analysis of 
the results, as are the errors on the displayed fits. 



aX 


FIG. 2. Phase shift as a function of A for m,, ~ 470 MeV. 


IV. RESULTS 


Fig. 0 shows an example plot of the ratio in Eq. ( [To| 
and the corresponding effective phase defined in Eq. (13) 
for A = 0.03. We observe a clear plateau in the effective 
phase for the illustrated fitting region, and correspond¬ 
ing linear behaviour in the ratio. As an aside, the fit 
indicates a clearly non-zero value for the t = 0 intercept, 
confirming that there is a small but statistically signifi¬ 
cant imaginary shift in the wavefunction overlap factors 
(given in Eq. 0 ) for this value of A. 

Repeating this procedure for each value of A and ex¬ 
tracting the phase shifts, we are able to calculate the 
linear shift with respect to A, illustrated in Fig.[^ From 
Eq. (§ we know that this shift is directly proportional to 
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FIG. 3. Effective phase plots for A = —0.025, 0.05 respec¬ 
tively at niTr ~ 310 MeV. The results in the second plot have 
greater statistics by a factor of 4. Note that the sign of the 
fitted value is highly dependent on the fit window. 

the disconnected contribution to AS. Since there is no 
phase shift in the zero-field limit, we have used a single¬ 
parameter slope fit to extract the linear shift. This anal¬ 
ysis is repeated at the lighter pion mass. Table |T] includes 
the calculated phase shift for each value of A on the en¬ 
sembles generated, and results of the described analyses 
are summarised in Table [iTl Using the methods outlined 
in [23], we have also calculated the individual connected 
contributions to AS on these ensembles, and hence are 
able to calculate the total (connected and disconnected) 
value of AS. 

At the lighter mass, we find a result consistent with 
zero for AEdisc.- This unusual result may be the result 
of a couple of different factors. The A values chosen may 
simply be too small, and the phase shift is not able to 
rise above the background correlator noise. Fig. show 
effective phase plots for the two values of A realised at 
this lighter quark mass, and show that there is no clearly 
identifiable plateau at these statistics. 

Alternatively, there may be a sign change in either the 
light or strange contribution to AS at some mass be¬ 
tween ttItt = 310 — 470 MeV. This is unlikely, however, 
as previous results at similar masses have shown a sig¬ 
nificant As contribution, which would require the light 
quark contribution to have a strong quark mass depen¬ 
dence. 

From |26j we have both non-singlet and singlet renor¬ 
malisation factors for the axial current at the SU(3) sym¬ 
metric point, 

^TNS = 0.8458(8), (14) 

= 0.8662(34). (15) 

Further calculations at additional quark masses are re¬ 
quired to perform a chiral extrapolation of these quanti¬ 
ties, however the pion mass dependence of these factors is 


expected to be mild based on the non-singlet calculation 
of Ref. mi. 

To obtain the renormalised total spin contribution we 
use the singlet renormalisation factor: 

AE^= zffAS^'"**'. ( 16 ) 

For the purely disconnected quantity, we include the mix¬ 
ing of the connected and disconnected contributions un¬ 
der renormalisation: 

A V'MS /vMS A I f '7MS /7MS A A /'T7\ 

^^disc. — ^A,S^^disc. + \^A,S ~ ^A,NS J ^^conn. • U 0 

Using the renormalisation factors from the SU(3) sym¬ 
metric point, we quote our MS results in the final two 
columns of Table ttn 

Finally, since at the SU(3) symmetric point all quarks 
contribute exactly the same amount to AE^®^ , then at 
this point we can determine As 

As^(to^ = 465MeV) = ^AS^^ = -0.018(6) . (18) 

O 

V. CONCLUDING REMARKS 

Culminating in the results of Table |ll| we have shown 
how the Feynman-Hellmann theorem may be applied to 
perform full calculations of hadronic matrix elements. 

Extensions of these particular calculations include 
higher-statistics simulations, particularly at the lighter 
pion mass, and the generation of ensembles at additional 
pion masses to identify the quark mass dependence of 
AS. Further analysis of the existing data should allow 
for the extraction of disconnected quark spin contribu¬ 
tions for the other octet baryons and the vector mesons. 

The FH technique demonstrated here can be easily 
adapted to a variety of other disconnected quantities, 
such as the gluonic contribution to angular momentum, 
which would otherwise be rather challenging using con¬ 
ventional approaches. 
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(ki, Ks) 


A Y^latt. 
^^disc. 


A j,MS(2 GeV^) 

■^disc. 

^jiMS{2 GeV^) 

(0.120900,0.120900) 

(0.121095,0.120512) 

1.001(7) -0.310(5) 
1.004(10) -0.319(6) 

-0.079(21) 0.612(24) 
0.014(16) 0.699(25) 

-0.055(18) 

0.026(14) 

0.530(21) 

0.605(21) 


TABLE II. Table of results at each pion mass for the individual quark axial charges and the disconnected and full (connected 
plus disconnected) contribution to the total quark spin. The quantities reported with the “latt.” superscript are unrenormalised. 
The final two columns report our renormalised results based upon Eqs. (141 through (17l. 
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